#export LD_PRELOAD=/usr/lib/libmpi.so python
from libatomism import *
from math import *
import array 
import matplotlib.pyplot as plt

def setrepro(kinop,N,Emin):

   gcoors=kinop.getCoordinates()
   unit = gcoors.getUnits()
   system = kinop.getSystem()
   gcoorsg = msGeneralizedCoordinates.New(unit)
   q0 = msScalarVariable.New("Angstrom").set(1., 0., 20, 0.05, 0.05).setId("q0")
   q1 = msScalarVariable.New("Angstrom").set(1., 0., 20, 0.05, 0.05).setId("q1")
   gcoorsg.addVar(q0).addVar(q1)
   
   repro =msRepromodel.New("")
   repro.setCoordinates(gcoorsg)
  
   d=0.1
   with open('PESdns.dat') as f:
     for line in f: 		# read rest of lines
        array = [float(x) for x in line.split()]
        if(len(array)>1):
	   value = (array[len(array)-1] + 188.656188) * 627.51
	   array.remove(array[len(array)-1])
	   if( (value>Emin) and(array[0]<2.3)):
	      vec = VectorOfDouble()
              vec.extend( arg for arg in array)
              gcoors.setCoordinates(vec)
              kinop.setDynamicDofs()
              system.computeCartCoordinates()
              vec = VectorOfDouble()
              vec.extend([system.separation(0,2)**(-2),system.separation(1,3)**(-2)])
              repro.addData( vec, value )
              
   msLogger.setPriorityMin(Priority.INFO)
   msLogger.setVerboseLevel(1)
   repro.setBasisSet(N)
   msLogger.setPriorityMin(Priority.ERROR)
   repro.fit()

   sampler=repro.getSampler()
   sampler.begin()
   field=repro.getField()
   y=[]
   x=[]
   n=0
   sigma=0
   while(sampler.getNextPoint(gcoorsg)):
      x.append(field.evaluate())
      y.append(repro.evaluate())
      sigma+=abs(field.evaluate() - repro.evaluate())
      n=n+1
   print "sigma=",sigma/n
   plt.plot(x,y,"+")
   plt.show()
   return repro
     
